function [ycounter,stcounter,countobs,countmat]=lmj_counterfactual_sub(funcmod,parvec,solveopt,addsol,countszero,etamat)
%% 1. Solution of the model
%     FIELD contains structures used to set-up extended state-space 
[G,R,Cs,eu,SDX,Zs,field]=feval(funcmod,parvec,solveopt,addsol);
if ~isequal(eu,[1;1]);
    ycounter=[];
    stcounter=[];
    countobs=[];
    countmat=[];
    return;
end
% =========================================================================
%% 2. Set-up extended states space 
%     See buildState.m by AB for an alternative 
%% 2.1 Dimensions and time invariant matrices 
[nobs nz]=size(etamat); 
ns       =size(G,1); 
% nzHF: Number of observables in High frequency
nzHF = field.nzHF;
if size(Zs,1)~=nzHF;error('Rows in Z must match number of variables observed in High frequency'); end 
% nzLF: Number of observables in Low  frequency 
nzLF = field.nzLF;
nmu  = size(SDX,2);
Z  = [Zs zeros(nzHF,nzLF); zeros(nzLF,ns) eye(nzLF)];
% Selection matrix 
A   = field.AA;
C   = [Cs;A*Cs]; 
% Number of periods of High frequency per 1 in Low frequency 
NUMIT= field.NUMIT; 
AggregationType = field.flag_tagsum;  % type of temporal aggregation
if AggregationType == 1
    sigma =[0 ones(1,NUMIT-1)];
else
    sigma =(1:NUMIT);
end
casevec =kron( ones(nobs/NUMIT,1) , (1:NUMIT)' );
%% 2.2 Page matrices with pages corresponding to quarters {1,2,3,4} 
Tpag  = zeros(ns+nzLF,ns+nzLF,NUMIT);
Rpag  = zeros(ns+nzLF,nmu    ,NUMIT);
for ii=1:NUMIT 
    if AggregationType == 1
        Tpag(:,:,ii) = [G zeros(ns,nzLF);A*G eye(nzLF)*sigma(ii)];
        Rpag(:,:,ii) = [R; A*R];
    else
        Tpag(:,:,ii) = [G zeros(ns,nzLF); (1/sigma(ii))*A*G eye(nzLF)*(sigma(ii)-1)/sigma(ii)];
        Rpag(:,:,ii) = [R; (1/sigma(ii))*A*R];
    end 
end 

%% 7.B. Decomposition of states and observables 
nstilda=ns+nzLF; 
nx=size(SDX,1); 
disp('Begin Counterfactual')
countmat=zeros(nobs,nstilda,nx);
countobs=zeros(nobs,size(Z,1),nx);

nz=size(Z,1); 
% Begin Counterfactual
for ii=1:nx
    tempeta=zeros(nx,nobs);
    tempeta(ii,:)=etamat(:,ii)';
    % ==================================================
    % Loop follows KFILTER_FULL.m
    ysim=zeros(nz,nobs);
    stsim=zeros(nstilda,nobs);
    stsim(:,1)=countszero(:,ii);
    ysim(:,1)=Z*stsim(:,1);
    for jj=2:nobs;
        stsim(:,jj)=Tpag(:,:, casevec(jj) )*stsim(:,jj-1) + Rpag(:,:,casevec(jj) )*(tempeta(:,jj));
        ysim(:,jj) =Z*stsim(:,jj);
    end
    % =====================================================
    countmat(:,:,ii)=stsim(:,:)';
    countobs(:,:,ii)=ysim';
end

ycounter=zeros(nz,nobs);
stcounter=zeros(nstilda,nobs);
stcounter(:,1)=sum(countszero,2);
ycounter(:,1)=Z*stcounter(:,1);
for jj=2:nobs;
    stcounter(:,jj)=Tpag(:,:, casevec(jj) )*stcounter(:,jj-1) + Rpag(:,:,casevec(jj) )*(etamat(jj,:)');
    ycounter(:,jj) =Z*stcounter(:,jj);
end
ycounter=ycounter';
stcounter=stcounter'; 

maxdif=comparemat( squeeze( sum(countobs,3) ),ycounter); 
if maxdif > 1e-5; 
    error('Cannot decompose observables across shocks') 
end 
dispaj('Maximum discrepancy Y and sum Y per shock=',maxdif); 
maxdif=comparemat( squeeze( sum(countmat,3) ),stcounter); 
if maxdif > 1e-5; 
    error('Cannot decompose states across shocks') 
end 
dispaj('Maximum discrepancy S_T and sum S_T per shock=',maxdif); 
clear stsim ysim; 
